Mapping guides: -More general intro: https://www.joshuastevens.net/cartography/make-a-bivariate-choropleth-map/ -R tutorial: https://cran.r-project.org/web/packages/biscale/vignettes/biscale.html

Simple Univariate Chloropleth Maps using Leaflet

Load libraries and dataset

library(tidyverse)
library(leaflet)
library(sf)
## Linking to GEOS 3.8.1, GDAL 3.1.1, PROJ 6.3.1
library(rgdal)
## Loading required package: sp
## rgdal: version: 1.5-18, (SVN revision 1082)
## Geospatial Data Abstraction Library extensions to R successfully loaded
## Loaded GDAL runtime: GDAL 3.1.1, released 2020/06/22
## Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/gdal
## GDAL binary built with GEOS: TRUE 
## Loaded PROJ runtime: Rel. 6.3.1, February 10th, 2020, [PJ_VERSION: 631]
## Path to PROJ shared files: /Library/Frameworks/R.framework/Versions/4.0/Resources/library/rgdal/proj
## Linking to sp version:1.4-4
## To mute warnings of possible GDAL/OSR exportToProj4() degradation,
## use options("rgdal_show_exportToProj4_warnings"="none") before loading rgdal.
# read in project data
merge_data = read_csv("data/merge_data.csv") %>% 
  rename(STUSPS = state_id) %>% # to match shapefile col name below
  mutate_if(is.numeric, round, digits = 2) # round to 2 decimal places
## 
## ── Column specification ────────────────────────────────────────────────────────
## cols(
##   .default = col_double(),
##   state_id = col_character(),
##   hostility = col_character()
## )
## ℹ Use `spec()` for the full column specifications.
merge_data
## # A tibble: 50 x 27
##    STUSPS counties_no_pro… percent_abortion percent_birth percent_medicaid
##    <chr>             <dbl>            <dbl>         <dbl>            <dbl>
##  1 AK                   79               11            73               24
##  2 AL                   91               10            74               18
##  3 AR                   97                9            75               28
##  4 AZ                   80               11            73               24
##  5 CA                   24               20            65               27
##  6 CO                   77               13            71               20
##  7 CT                   13               23            62               23
##  8 DE                   33               19            66               19
##  9 FL                   67               21            64               17
## 10 GA                   94               17            68               13
## # … with 40 more rows, and 22 more variables: percent_women_no_provider <dbl>,
## #   percent_uninsured <dbl>, percent_bc_18_49 <dbl>,
## #   abortion_rate_15_17_state <dbl>, abortion_rate_15_19_state <dbl>,
## #   abortion_rate_15_44_state <dbl>, abortion_rate_18_19_state <dbl>,
## #   birthrate_15_17_state <dbl>, birthrate_15_19_state <dbl>,
## #   birthrate_18_19_state <dbl>, need_bc_hisp_20_44 <dbl>,
## #   need_bc_hisp_younger_20 <dbl>, need_bc_13_44 <dbl>,
## #   need_bc_black_20_44 <dbl>, need_bc_black_under_20 <dbl>,
## #   need_bc_white_20_44 <dbl>, need_bc_white_under_20 <dbl>,
## #   need_bc_under_20 <dbl>, total_expend_abortion <dbl>,
## #   public_expenditures <dbl>, expenditure_rate <dbl>, hostility <chr>
# import state shapefile from data folder
all_states <- readOGR("data/cb_2018_us_state_500k/cb_2018_us_state_500k.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/rht2112/Documents/My Documents/Grad School/Coursework/Semester 3/Data Science I/Homeworks and Projects/p8105_2020_final/data/cb_2018_us_state_500k/cb_2018_us_state_500k.shp", layer: "cb_2018_us_state_500k"
## with 56 features
## It has 9 fields
## Integer64 fields read as strings:  ALAND AWATER
summary(all_states)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min       max
## x -179.1489 179.77847
## y  -14.5487  71.36516
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=NAD83 +no_defs]
## Data attributes:
##    STATEFP            STATENS            AFFGEOID            GEOID          
##  Length:56          Length:56          Length:56          Length:56         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##     STUSPS              NAME               LSAD              ALAND          
##  Length:56          Length:56          Length:56          Length:56         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##     AWATER         
##  Length:56         
##  Class :character  
##  Mode  :character

Join project data with sf file

data_sf <- merge(all_states, merge_data, all.x = F) 
summary(data_sf)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##          min       max
## x -179.14891 179.77847
## y   18.91036  71.36516
## Is projected: FALSE 
## proj4string : [+proj=longlat +datum=NAD83 +no_defs]
## Data attributes:
##     STUSPS            STATEFP            STATENS            AFFGEOID        
##  Length:50          Length:50          Length:50          Length:50         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     GEOID               NAME               LSAD              ALAND          
##  Length:50          Length:50          Length:50          Length:50         
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##                                                                             
##     AWATER          counties_no_provider percent_abortion percent_birth  
##  Length:50          Min.   :13.00        Min.   : 5.00    Min.   :57.00  
##  Class :character   1st Qu.:75.50        1st Qu.: 9.00    1st Qu.:68.00  
##  Mode  :character   Median :91.00        Median :11.50    Median :72.50  
##                     Mean   :78.24        Mean   :13.36    Mean   :71.12  
##                     3rd Qu.:95.00        3rd Qu.:17.00    3rd Qu.:75.00  
##                     Max.   :99.00        Max.   :29.00    Max.   :79.00  
##                                                                          
##  percent_medicaid percent_women_no_provider percent_uninsured percent_bc_18_49
##  Min.   : 9.00    Min.   : 0.00             Min.   : 3.00     Min.   :61.80   
##  1st Qu.:15.50    1st Qu.:18.25             1st Qu.: 7.00     1st Qu.:67.90   
##  Median :19.50    Median :42.50             Median :11.00     Median :71.10   
##  Mean   :20.24    Mean   :43.16             Mean   :10.82     Mean   :70.44   
##  3rd Qu.:24.75    3rd Qu.:62.75             3rd Qu.:13.00     3rd Qu.:73.45   
##  Max.   :38.00    Max.   :96.00             Max.   :24.00     Max.   :78.00   
##                                                               NA's   :12      
##  abortion_rate_15_17_state abortion_rate_15_19_state abortion_rate_15_44_state
##  Min.   :1.000             Min.   : 2.600            Min.   : 4.600           
##  1st Qu.:2.200             1st Qu.: 4.600            1st Qu.: 7.675           
##  Median :3.050             Median : 6.000            Median : 9.900           
##  Mean   :3.468             Mean   : 6.684            Mean   :11.238           
##  3rd Qu.:4.450             3rd Qu.: 8.200            3rd Qu.:14.100           
##  Max.   :9.600             Max.   :17.200            Max.   :28.200           
##                                                                               
##  abortion_rate_18_19_state birthrate_15_17_state birthrate_15_19_state
##  Min.   : 4.80             Min.   : 3.500        Min.   : 8.50        
##  1st Qu.: 8.20             1st Qu.: 6.400        1st Qu.:15.93        
##  Median :10.00             Median : 8.250        Median :19.75        
##  Mean   :11.49             Mean   : 8.648        Mean   :20.72        
##  3rd Qu.:13.25             3rd Qu.:10.175        3rd Qu.:24.80        
##  Max.   :28.40             Max.   :15.000        Max.   :34.70        
##                                                                       
##  birthrate_18_19_state need_bc_hisp_20_44 need_bc_hisp_younger_20
##  Min.   :13.80         Min.   :    680    Min.   :   270         
##  1st Qu.:29.75         1st Qu.:  11792    1st Qu.:  2745         
##  Median :37.85         Median :  28845    Median :  7275         
##  Mean   :38.99         Mean   :  82409    Mean   : 19617         
##  3rd Qu.:48.75         3rd Qu.:  66735    3rd Qu.: 16162         
##  Max.   :66.50         Max.   :1102300    Max.   :263760         
##                                                                  
##  need_bc_13_44     need_bc_black_20_44 need_bc_black_under_20
##  Min.   :  34960   Min.   :   160      Min.   :   60         
##  1st Qu.: 119260   1st Qu.:  4402      1st Qu.: 1148         
##  Median : 303910   Median : 26175      Median : 6650         
##  Mean   : 411957   Mean   : 59597      Mean   :14382         
##  3rd Qu.: 463600   3rd Qu.:101400      3rd Qu.:24575         
##  Max.   :2526010   Max.   :264700      Max.   :61270         
##                                                              
##  need_bc_white_20_44 need_bc_white_under_20 need_bc_under_20
##  Min.   : 11630      Min.   :  2310         Min.   :  8220  
##  1st Qu.: 52142      1st Qu.: 18342         1st Qu.: 26780  
##  Median :121200      Median : 41105         Median : 68105  
##  Mean   :148834      Mean   : 52164         Mean   : 92522  
##  3rd Qu.:213080      3rd Qu.: 71735         3rd Qu.:103850  
##  Max.   :468270      Max.   :161870         Max.   :527440  
##                                                             
##  total_expend_abortion public_expenditures expenditure_rate  hostility        
##  Min.   :    0         Min.   :   762      Min.   :0.0100   Length:50         
##  1st Qu.:    0         1st Qu.:  8054      1st Qu.:0.0500   Class :character  
##  Median :   76         Median : 24971      Median :0.0900   Mode  :character  
##  Mean   : 1931         Mean   : 41698      Mean   :0.0924                     
##  3rd Qu.:  453         3rd Qu.: 56273      3rd Qu.:0.1200                     
##  Max.   :32613         Max.   :454706      Max.   :0.2700                     
##  NA's   :13

Abortion rates map

# make color palette
abortion_rate_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$percent_abortion,
  reverse = TRUE)

# base map
base_map <- leaflet(data_sf) %>%
  addProviderTiles("CartoDB.PositronNoLabels") 

# add in custom labels for percent abortion
labels_abortion_rate <- sprintf(
  "<strong>%s</strong><br/>%g&#37;",
  data_sf$STUSPS, data_sf$percent_abortion
) %>% lapply(htmltools::HTML)

# full map
abortion_map = base_map %>% 
  addPolygons(
  fillColor = ~abortion_rate_pal(data_sf$percent_abortion),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_abortion_rate,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = abortion_rate_pal, values = ~percent_abortion, opacity = 0.7, title = 'Abortion Rates by State <br> (% of pregnancies ending in abortions)',
  position = "bottomright")
abortion_map

Public contraceptive expenditure rate map

# make color palette
bc_fund_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$expenditure_rate,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_bc_fund <- sprintf(
  "<strong>%s</strong><br/>&#36;%g/woman",
  data_sf$STUSPS, data_sf$expenditure_rate
) %>% lapply(htmltools::HTML)

# full map
bc_fund_map = base_map %>% 
  addPolygons(
  fillColor = ~bc_fund_pal(data_sf$expenditure_rate),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_bc_fund,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = bc_fund_pal, values = ~expenditure_rate, opacity = 0.7, title = 'Public Contraceptive Expenditure Rate by State  <br> ($/woman in likely need of publicly-funded services)',
  position = "bottomright")

bc_fund_map

Abortion access map

# make color palette
abortion_access_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$percent_women_no_provider,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_abortion_access <- sprintf(
  "<strong>%s</strong><br/>%g&#37;",
  data_sf$STUSPS, data_sf$percent_women_no_provider
) %>% lapply(htmltools::HTML)

# full map
abortion_access_map = base_map %>% 
  addPolygons(
  fillColor = ~abortion_access_pal(data_sf$percent_women_no_provider),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_abortion_access,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = abortion_access_pal, values = ~percent_women_no_provider, opacity = 0.7, title = 'Percent of women living in a county <br> without access to abortion clinic by State',
  position = "bottomright")

abortion_access_map

Teen birth rates map

# make color palette
teen_births_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$birthrate_15_19_state,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_teen_births <- sprintf(
  "<strong>%s</strong><br/> %g births per 1000 women",
  data_sf$STUSPS, data_sf$birthrate_15_19_state
) %>% lapply(htmltools::HTML)


# full map
teen_births_map = base_map %>% 
  addPolygons(
  fillColor = ~teen_births_pal(data_sf$birthrate_15_19_state),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_teen_births,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = teen_births_pal, values = ~birthrate_15_19_state, opacity = 0.7, title = 'No. of births per 1000 women aged 15-19',
  position = "bottomright")

teen_births_map

Number of Women under 20 in need of accessing publicly-funded contraceptives

# make color palette
need_bc_teen_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$need_bc_under_20,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_need_bc_teen <- sprintf(
  "<strong>%s</strong><br/>%g",
  data_sf$STUSPS, data_sf$need_bc_under_20
) %>% lapply(htmltools::HTML)

# full map
need_bc_teen_map = base_map %>% 
  addPolygons(
  fillColor = ~need_bc_teen_pal(data_sf$need_bc_under_20),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_need_bc_teen,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = need_bc_teen_pal, values = ~need_bc_under_20, opacity = 0.7, title = 'No. of women who likely need public support for contraceptive services and supplies younger than 20',
  position = "bottomright")

need_bc_teen_map

Total Number of Women in need of accessing publicly-funded contraceptives

# make color palette
need_bc_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$need_bc_13_44,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_need_bc <- sprintf(
  "<strong>%s</strong><br/>%g",
  data_sf$STUSPS, data_sf$need_bc_13_44
) %>% lapply(htmltools::HTML)

# full map
need_bc_map = base_map %>% 
  addPolygons(
  fillColor = ~need_bc_pal(data_sf$need_bc_13_44),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_need_bc,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = need_bc_pal, values = ~need_bc_13_44, opacity = 0.7, title = 'Total No. of women who likely need public support for contraceptive services and supplies',
  position = "bottomright")

need_bc_map

State expenditure rate

# make color palette
expenditure_rate_pal <- colorNumeric(
  palette = "inferno",
  domain = data_sf$expenditure_rate,
  reverse = TRUE)

# add in custom labels for public contraceptive expenditure
labels_expenditure_rate <- sprintf(
  "<strong>%s</strong><br/>%g",
  data_sf$STUSPS, data_sf$expenditure_rate
) %>% lapply(htmltools::HTML)

# full map
expenditure_rate_map = base_map %>% 
  addPolygons(
  fillColor = ~expenditure_rate_pal(data_sf$expenditure_rate),
  weight = 2,
  opacity = 1,
  color = "white",
  dashArray = "3",
  fillOpacity = 0.7,
  highlight = highlightOptions(
    weight = 5,
    color = "#666",
    dashArray = "",
    fillOpacity = 0.7,
    bringToFront = TRUE),
  label = labels_expenditure_rate,
  labelOptions = labelOptions(
    style = list("font-weight" = "normal", padding = "3px 8px"),
    textsize = "15px",
    direction = "auto")) %>% 
  addLegend(pal = expenditure_rate_pal, values = ~expenditure_rate, opacity = 0.7, title = 'Ratio of Total reported public expenditures for family planning client services in 1000s of dollars to number of women in likely need of public support for contraceptive services and supplies',
  position = "bottomright")

expenditure_rate_map

Hostility map